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m ; Abstract 

S-^ ■ In this work, we carefully study the energy eigen-values and splitting 

Q"^ , of heavy quarkonia as there exist 1/r^ and (5^(r) singular terms in the 

0^ I potential which make a direct numerical solution of the Schrodinger 

^S^ • equation impossible. We compare the results obtained in terms of per- 

^^. turbation and variational methods with various treatments. 

Oh! 



X 



It is well known that the interactions between heavy quarks (or heavy quark-antiquark) can 
be reasonably described by a combination of the Coulomb-type potential which is induced by 
the single-gluon exchange and the confinement potential. That form corresponds to the zero-th 
order potential. 

As the first order corrections, the Breit-Fermi potential can be derived directly from the 
elastic scattering amplitude of the quark and antiquark, where the non-relativistic approximation 
is employed jj|]. 

In the derivation, a non-relativistic expansion with respect to (|p|/Mq)" (n > 1) powers is 

carried out and usually, one only keeps terms up to \p\'^/Mq and ignore the higher orders. A 

Fourier transformation would produce 5'^{f) and 1/r^ terms in the configuration space. The 

extra potential term which is related to quark spins can be written in a general form as |^] 
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where S* = si + 5*2, r is a unit spatial vector, L is the relative orbital angular momentum of the 
quark-antiquark system. 

As noted, even though the Coulomb potential 1/r is also singular at r ^ 0, it is benign for 
solving the second-order differential equation. In fact, the Coulomb potential corresponds to 
the zero-th order of the non-relativistic expansion of the one-gluon-exchange induced potential, 
and its singular behavior is fully compensated by the small measure at vicinity of r ^ 0. By 
contraries, the extra 5^(r) and 1/r^ come from higher order (at least the first order) of the 
expansion. They are more singular and make the equation unsolvable. Here the unsolvability 
means that when one tries to numerically solve the Schrodinger equation with the potential 
including 1/r'^ and 6^{r) in expression (|^), he would confront infinity, namely the solution blows 
up due to the singular behavior of the potential at vicinity of r — > 0, i.e., no physical solution can 
be obtained. Landau pointed out that as the potential is more singular than 1/r^, the equation 
is unsolvable [0], obviously 1/r^ and 5'^{f) are more singular than 1/r^. 

Analysis tells us that such high order singularities are not physical, but brought up by 
artificial truncation of the non-relativistic expansion, namely one only keeps to \p\'^/Mq (or 
even higher, but finite orders). 

Gupta et al. suggested to take an alternative expansion with respect to \p\/E instead of 
\p\/Mq 0]. Obviously the divergence for |p| ^ 00 (or r ^ 0) is avoided, because \p\/E — > 1 as 
IpI -^ cxD. However, E is an operator itself, while it exists at denominator of the expression, the 
equation is not solvable either unless one takes some extra approximations Q. 

As discussed above, the singularity is artificial because of the non-relativistic approximation, 
so that in principle, if one deals with the problem in a reasonable scheme, for example, the 
relativistic B-S equation, such singularity would not appear at all. However, the potential form 
which is non-relativistic, is most useful and convenient. Many good zero-th order results have 
been achieved H, we have all reasons to believe that the non-relativistic Schrodinger equation 
properly describes the physics. When the first order correction related to quark spins are 



considered, which determines the fine-structure and hyperfine structure of the quarkonium, the 
singularity problem emerges, and it turns the whole equation unsolvable. Therefore, as long as we 
are working with the non-relativistic Schrodinger system, we are forced to face this unphysical 
singularity problem. Similar singularity problem was discussed in previous literatures. For 
example, Frank, Land and Spector analyzed the singular potential [^, however, they mainly 
concerned the scattering states instead of the bound states. In this work, we only focus on 
the bound state problems. Our goal is to correctly solve the problem and obtain the physical 
quantities, even though there is the singularity. 

So now the question is can we obtain the eigen energies and wavefunctions which are suf- 
ficiently close to physical reality? Of course, the closeness or validity degree can be tested by 
comparing the obtained results about the splitting, for example, with experimental data. 

The effect of the 6^{f) term on the s-wave has been carefully studied by some authors [^, so 
readers are recommended to the original work for details. 

In this work, we only concentrate ourselves at 1/r^ singular term which does not affect the 
s-wave (see below). 

In this paper, we will obtain the eigen energies in both perturbative and variational meth- 
ods and then compare the results, we also discuss some alternative ways to modify the potential 
form for avoiding the singularity which have been employed in literatures. Our numerical results 
show that there exist certain regions for those possible parameters in the approaches where all 
the results are consistent to each other with sufficient accuracy. 

(i). Perturbative method. 

The zero-th order potential does not include the singular parts of (|l|), so that we can easily 
solve the corresponding Schrodinger equation and obtain the zero-th order eigen energies and 
wavefunctions. Then we are going on studying the first order equation which involves the singular 
terms and obtain corresponding energies in terms of perturbative method. 

Let us investigate the singular terms related to 1/r^ in eq.(||), we find that L ■ S = ^[{L + 
S)^ — L^ — S'^], so for the S-state (/ = 0), its eigenvalue is zero, moreover, the tensor-coupling 



term in (|^) is proportional to an operator 

Vtensor = [3(si • r)(s2 • f) - (si • S2)], 

which is sandwiched between two states which are solutions of the zero-th order equation. If 
both sides are S-states, then 

dn < loolV^tensorl^^OO >= 0, 

where ri is the solid angle. The lowest non-zero value comes from P-state, < / = 1| VjensorK — 
1 >7^ 0, and the operator can also result in a small mixing between S-state and D-state, i.e. 
< I = 0\Vtensor\^ = 2 >7^ 0. We are not going to discuss the details of the mixing in this work. 

Therefore, the singular 1/r^ terms do not bother the S-states, but only affect the higher 
orbital-angular momentum states {I > 1), this enables us to avoid singularities in practical 
calculations (see below). With the first order correction the Schrodinger equation for / = 1,5"= 1 
and J = 2, 1, {S = si + 5*2) is recast as 

where we only keep the L ■ S term (dropped out the tensor coupling for our purpose of this 
study) and take mi = m2 = m = Mq, the reduced mass /i = m/2. 

In this study, our purpose is to investigate the effects of the singular potential, so that in the 
Schrodinger equation we deliberately omit the confinements piece, such Kr which is well behaved 
at r — > for convenience. This treatment does not change the results of the study, but greatly 
simplifies the calculation. 

(ii) The practical calculation in terms of perturbation. 

For heavy quarkonia cc, bb, the coefficients of the 1/r^ is much smaller than that of 1/r 
which has a milder behavior at r ^ 0, so one can employ the perturbation method to calculate 
the splitting. Then the zero-th order equation is 

^v'-^T^^'t^ + a - ^)R^'Hr) = E(^)R(^Hr). (3) 
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For the Ip state, there is an analytical solution as 



and the Bohr radius is 
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If we take the conventional parameters rric = 1.84 GeV, as = 0.39, then 

a = 2.09030 GeV"\ E^^'^ = -0.031096 GeV. 

We can evaluate the contribution of the 1/r^ term by means of perturbation. The perturbative 
hamiltonian is 

Using the zero-th order wavefunction (^), we have 

Ip = <P\5H\P>o, j Rf{r)^(^)^Rf{rydr 

It is worth of noting, r exists at numerator of the w~^'{r) expression in (^), so that the 
integration converges. 

Numerically, we obtain 

Ip = 0.00105105 GeV. (9) 

For 2^P2, 2^Pi and 2'^Pq states, the energy shifts are Ip, —Ip and —2Ip respectively. These 
values are only 3~ 6% of E^^' , therefore the perturbation is legitimate and one can be convinced 
that the results are reliable. 

(iii) Variational method. 

The Schrodinger equation is not solvable if there is a potential term —r~^ with s > 2 0] 
or (5^(r), it is because the corresponding eigen-energy tends to -co and average radius turns 
to 0. "Particle falls to the point r = 0". Let us focus on the 1/r^ singularity. Definitely, the 



— oo solution is not physical, but a virtual solution of the equation with the singular potential. 
The key point is how to systematically get rid of the virtual unphysical solution and achieve the 
physical one. 

The hamiltonian is 
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where we take J = 1 as an example. If we employ the P-state wavefunction (^ as a trial 
function, the average value is 

where a is a variational parameter in the expression. Later we need to differentiate the expression 
with respect to a and obtain the minimum. 

For a cc bound state, the relation between < hp > and a shows that as a ^ 0, < hp >^ — oo, 
there is no lower bound. The corresponding average radius is 

< r >= 5a, (12) 

so that as a ^ 0, < r >^ 0. It indicates that if one took the minimal value — oo as the binding 
energy, the radius turns out to be infinitesimal, it cannot exist in principle. It is exactly the 
Landau's conclusion. We can also prove that the conclusion is independent of the adopted form 
of the trial functions. This result indicates that the singular behavior of the potential near r — > 
is transferred to < /ip > at a ^ 0. 

Strictly speaking, the equation is of no solution in this case, does it mean that the Schrodinger 
equation cannot describe our physics picture? As pointed out in previous section, this singu- 
larity is caused by the ill behavior of the non-relativistic approximation near r = 0. One can 
conjecture that the given Schrodinger equation can still correctly describe the physics even with 
the unphysical singular terms, but we have to solve it in a proper way. Namely since the singu- 
larity at r — > is artificial and unphysical, the boundless energy corresponding to a = does 
not reflect physical reality, i.e. the infinite < /ip > at a = is not a physical solution and must 
be eliminated. 



From another angle, as o ^ 0, the trial function (^) has a limit as 

lim -^e-''/" = 0. (13) 

The trial function tends to zero, the solution becomes trivial and does not make any sense at 
ah. 

In fact, if one differentiate < hp > with respect to a, besides the minimal value — oo at the 
boundary a = 0, < hp > has a local minimum at oi and it is 

a.= '^+'^"^-"^°---="»l.!,7832GeV-. (14) 

32as 

The corresponding average energy value is 

Eyar = -0.0322362 GeV. (15) 

When we compare this value Eyar with the result obtained in terms of perturbation 

Eper = E^^^ -Ip = -0.0321470 GeV, (16) 

the deviation between E^ar and -Eper is only about 10~^. 

For a real meson, < r > cannot be zero, so the variational parameter a must be non-zero. 
With this constraint, the boundless energy which corresponds to a — > 0, is ruled out and we can 
enjoy the reasonable solution determined by the local minimum. Actually, this local minimum 
should correspond to the real physical value. 

Even though we know the negative infinite minimal value of E is not physical and should be 
abandoned, in practical calculations on computer, it is still bothering, therefore some plausible 
ways which may eliminate the singularity from the concerned problem would be welcome. 

(iv) Variational method and a "smear" form factor. 

One may introduce a form factor to "smear" the singularity Q. The most widely adopted 
form factor is the error function er/(Ar). It turns the spatial part of the spin-orbit coupling to 

Vr^=''-^erfiXr)[''h (17) 
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For any finite A— value, Vp is no longer singular as 
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which is finite as long as A is finite. In the approach, the hamiltonian turns into a form 
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Because it has a nice behavior at r ^ 0, in < /ip > there is no more the troublesome negative in- 
finity at a — > 0, and one can safely handle it with the variational method. For A = 5, 10, 15 GeV, 
we obtain the minimum of < hp > which is also the minimal value for the whole range of 
< a < cxD, and tabulate the results in table. 1. 



Table 1. 



A (GeV) 


5 


10 


15 


a (GeV-^) 


1.97967 


1.97868 


1.97849 


E (GeV) 


-0.0322276 


-0.322339 


-0.0322352 



For A = 15 GeV, the result is exactly the same as E^ar achieved in last section. But as 
A > 15 GeV, < hp > has a boundless negative minimal value for a — > again. It means that 
as A in the form factor is greater than 15 GeV, the function cannot smear out the singularity 
behavior at r ^ 0. The meaning of introducing the "smear" form factor is that once we intro- 
duce a " smear form factor whose parameter resides within a certain range" , we do not need to 
worry about the boundless minimal value at all, namely then the local minimum is the mini- 
mal value and finite, so that we can safely use the variational method to obtain energies. But 
beyond the parameter range, the boundless minimal value appears again, then we need to look 
for the local minimum as the real solution, exactly in analog to the treatment of previous section. 



(v) The variational method and the " regularization" . 

(a) Another approach to handle the singularity is the so-called "regularization" [p. In 
this approach, one can "regularize" pr to ,j.}^\n (c > 0), which obviously makes the potential 



converge at r ^ 0. Thus Vp is replaced by 



and V, 
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Vp as c ^ 0. We tabulate the results corresponding to various c— values in Table. 2. 



Table 2 



c (GeV-1) 


0.1 


0.05 


0.035 


a (GeV-^) 


1.99094 


1.98525 


1.98333 


E (GeV) 


-0.0321330 


-0.0321808 


-0.0321963 



It is noted that the smaller c is, the closer to Eyar and values in Table. 1 the i?— values 
obtained by the " regularization" scheme are. But as c < 0.03 GeV~^, the singular behavior of 

(r) 

the potential shows up again, namely then < hp > would have a boundless negative minimal 
value at a ^ 0. 

In fact, in sections (iv) and (v), two different approaches are introduced to remedy the 
singular behavior of the potential at r = 0. Besides the variational parameter a, there is an 
extra parameter A or c being introduced. We notice that if the concerned parameter falls in a 
certain range, the negative infinite minimal value oi < hp > which corresponds to an extreme 
variational parameter (usually a ^ 0) does not exist at all and a local minimum would give 
the physical solution. However, once the parameter is beyond the allowed range, such negative 
infinite < hp > would appear again, in this case, one has to enforce a constraint condition as in 
section (iii) for the pure variational method, i.e. prior rule out the negative unphysical infinity 
and keep the local minimum as reasonable solution. 

(b) As the simplest "regularization" scheme, one can replace Vp in eq.(pO|) by 

3 4a, 1 



VJ,''^ = 9{r - roh 



(21) 



' 2m'^ 3 r^ 

Obviously, it is an equivalent way to eq.(^) to get rid of the singularity at r ^ and tq can 
be seen as an arbitrarily introduced "cut-ofF' in the configuration space or a corresponding 
Aq ~ 1/ro in momentum space. For the zeroth order l^Pi state the correction of Vp can be 



calculated and the numerical results are shown in Table. 3. 



table 3. 



ro (GeV-^) 


0.20 


0.15 


0.10 


0.08 


0.068 


a (GeV-^) 


1.97926 


1.97886 


1.97857 


1.97848 


1.97844 


E (GeV) 


-0.0322303 


-0.0322328 


-0.0322347 


-0.0322352 


-0.0322355 



To make more sense, it would be interesting to re-calculate the energy corrections owing to 
the potential Vf^ for the l^Pi state in perturbation and compare 

'■°° , , -3 4a, 1 
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Now let us tabulate the results in Table 4 as 



(ii(0)) 



2 2i 



(22) 



(23) 



Table 4. 



ro (GeV-i) 


0.2 


0.15 


0.1 


0.08 


0.068 


E"-^ (GeV) 


-0.0010465 


-0.0010485 


-0.0010499 


-0.0010503 


-0.0010505 


5 (GeV) 


-4.51 X 10"" 


-2.58 X 10"'^ 


-1.17 X lO"*' 


-7.5 X 10-^ 


-5.4 X 10"^ 


R 


0.0043 


0.0024 


0.0011 


0.0007 


0.0005 



In the table, the ratio R is defined as 



R 



5 

E^' 



The integration over r from to tq denotes the part we ignore when taking the "regulariza- 
tion" (PI) with the step function 6{r — ro). The ratio R reflects the relative errors, and from 
Table 4, we can see the errors brought by the " regularization" are about a few thousandths 
at most. Therefore, according to the experimental accuracy, the scheme is very satisfactory. 
Indeed, our results restrict the probable error ranges for using the "regularization" or "smear" 
scheme to remedy the singularity. 
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(vi) With the Gupta's approach. 

In Gupta's approach 0], where an approximation s = p' +p being small is taken, p,p' are 
the three- momenta of the two quarks. Vp is replaced by 

where 

/i(r) = [1 - (1 + 2mr)e-2"'-]/mV2, (25) 

obviously 

lim/i(r) = 2. (26) 

Thus this modified Vp only possesses the singular behavior of 1/r, so that would not bring up 
any difficulties for solving the Schrodinger equation by variational method. With this approach, 
we employ the same trial function (Q) and variational parameter a to obtain energy as 

4^), = 1.98570 GeV-\ and E^^) = -0.0321867 GeV. (27) 

Instead, if we use the perturbation method, we obtain 

S^£^ = -0.0321085 GeV. (28) 

This is consistent with the values given in Table. 1 and 2, as the errors are negligibly small. 

(vii) Conclusion and discussion. 

Based on the above calculations, we can draw our conclusion as 

(a) As Landau pointed out that an ill-extrapolation of the non-relativistic expansion brings 
up an artificial singular r~^ which makes the Schrodinger equation unsolvable. Here the "un- 
solvable" means that it is impossible to obtain an exact solution of the equation, so that one is 
forced to look for a reasonable way to find the physical solution from the given hamiltonian with 
such singular terms. Our purpose is to seek for a solution which would describe the concerned 
physics (eigen-energies and eigen-wavefunctions) to a satisfactory degree. The success of the 
potential model for charmonia encourages this approach and stimulates further studies. 
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(b) In particular, the variational method without any special treatments for the singularity 
would result in an infinite negative energy which corresponds to < hp >= —oo, a = 0, < r >= 
0. Obviously it is not physical, because it manifests itself as an infinitesimal radius meson, i.e. 
a point-like particle which would cause singularity as well known. Moreover, as the variational 
parameter turns to zero, the trial function (^) is zero, so it becomes trivial and unphysical, so 
must be ruled out. 

The local minimum which corresponds to non-zero variational parameter and average radius 
indeed gives rise to reasonable solution and should be physical. Our results show a satisfactory 
consistency with that by other approaches, as expected. 

(c) The "smear" and "regularization" approaches all introduce new parameters whose phys- 
ical meaning is to decrease the measure near vicinity of r ^ and/or simply eliminate a small 
neighborhood of r = 0. Such schemes ("smear" or "regularization") with certain parameters 
partly compensate the artificial truncation of the non-relativistic expansion. Our results show 
that for certain parameter ranges, the obtained results are consistent with the perturbation. 
For future application, one needs to adopt reasonable values for the parameters, concretely, 
A = 15 GeV for the "smear" approach and c = 0.035 GeV~^ and vq = 0.068 for the "regulariza- 
tion". 

Beyond the parameter ranges, the singular behavior appears again. 

It may only concern the automatic calculation procedure on computer for obtaining the 
energies by variational method. Without carefully programming, the computation may overflow 
once it gets to the artificial unphysical a — > and refuses to look for the physical local minimum. 
Therefore when one writes his program, he has to be careful about the possible negative infinity 
solution. 

But if we take the suggested schemes with the concerned parameters residing in the allowed 
regions, we do not need to worry about the unphysical infinity at all while programming. 

(d) The wavefunction problem is not solved yet. One can use the perturbation method to 
obtain the first order wavefunction, but it needs to sum over contributions from all possible 
zeroth order states, so practical applications are restricted. With the variational method one 
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can obtain wavefunctions easily, but as we studied in another work Q, it is not a simple job. 
Even without the troublesome 1/r^ terms, it is difficult to obtain wavefunctions with satisfactory 
accuracy. Corresponding to the well-known potential forms, the Cornell, logarithmic etc. we 
have found that it is easy to obtain very accurate energy levels, but the errors for wavefunctions, 
in particular, for their values at origin are rather large unless we carefully choose the forms 
of the trial functions and use the multi-variational parameters. Moreover, the form of trial 
functions and number of variational parameters also depend on the form of potential. In |^ we 
only succeeded to select proper trial functions and number of parameters for each well-known 
potential forms. However, that is only for the zeroth order situation, so when we consider the 
first order correction in the potential, the situation is even worsened. 

So how to determine the wavefunctions, especially the values at origin, with existence of Xjr^ 
terms, is an open question and we will investigate it in our future works. 
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